Seismic and electromagnetic controlled-source interferometry in dissipative media

Seismic interferometry deals with the generation of new seismic responses by crosscorrelating existing ones. One of the main assumptions underlying most interferometry methods is that the medium is lossless. We develop an ‘interferometry-bydeconvolution’ approach which circumvents this assumption. The proposed method applies not only to seismic waves, but to any type of diffusion and/or wave field in a dissipative medium. This opens the way to applying interferometry to controlledsource electromagnetic (CSEM) data. Interferometry-by-deconvolution replaces the overburden by a homogeneous half space, thereby solving the shallow sea problem for CSEM applications. We demonstrate this at the hand of numerically modeled CSEM data. I N T R O D U C T I O N Seismic interferometry is the branch of science that deals with the creation of new seismic responses by crosscorrelating seismic observations at different receiver locations. Since its introduction around the turn of the century, the literature on seismic interferometry has grown spectacularly. Interferometric methods have been developed for random fields (Larose et al. 2006; Gerstoft et al. 2006; Draganov et al. 2007) as well as for controlled-source data (Schuster and Zhou 2006; Bakulin and Calvert 2006). The underlying theories range from diffusion theory for enclosures (Weaver and Lobkis 2001), stationary phase theory (Schuster et al. 2004; Snieder 2004) to reciprocity theory (Wapenaar et al. 2004; Weaver and Lobkis 2004; van Manen et al. 2005). All these theories have in common the underlying assumption that the medium is lossless and non-moving. The main reason for this assumption is that the wave equation in lossless non-moving media is invariant for time-reversal, which facilitates the derivation. Until 2005 it was commonly thought that time-reversal invariance was a necessary condition for interferometry, but recent research shows that this assumption can be relaxed. Slob, Draganov and Wapenaar (2006) analyzed the interferometric method for ground-penetrating radar data (GPR), in which losses play a prominent role. They showed that losses lead to amplitude errors as well as to the occurrence of spurious events. By choosing the recording locations in a specific way, the spurious events arrive before the first desired arrival and can thus be identified. Snieder (2006, 2007) followed a different approach. He showed that a volume distribution of uncorrelated noise sources, with source strengths proportional to the dissipation parameters of the medium, precisely compensates for the energy losses. As a consequence, the responses obtained by interferometry for this situation are free of spurious events and their amplitudes decay the way they should in a dissipative medium. This approach does not only hold for waves in dissipative media, but also for pure diffusion processes. Time-reversal invariance as well as source-receiver reciprocity break down in flowing or rotating media, but with some minor modifications interferometry also appears to work for these situations (Wapenaar 2006; Godin 2006; Ruigrok, Draganov and Wapenaar 2008). Recently we showed that interferometry, including its extensions for waves and diffusion in dissipative and/or moving media, can be represented in a unified form (Wapenaar et al. 2006; Snieder, Wapenaar and Wegler 2007). In turn, from this unified formulation it follows that the interferometric method can also be used for C © 2008 European Association of Geoscientists & Engineers 419 420 K. Wapenaar, E. Slob and R. Snieder more exotic applications like electroseismic prospecting and quantum mechanics. Interferometry in the strict sense makes use of crosscorrelations, but in the following we will extend the definition of interferometry so that it also includes crossconvolution and deconvolution methods. Slob, Draganov and Wapenaar (2007) introduce interferometry by crossconvolution and show that it is valid for arbitrary dissipative media. The crossconvolution method does not require a volume distribution of sources, but one restriction is that it only works for transient signals in specific configurations with receivers at opposite sides of the source array. The latter restriction does not apply to ‘interferometry-by-deconvolution’, which is the method discussed in this paper. I N T E R F E R O M E T RYB YD E C O N V O L U T I O N : 1 D V E R S I O N ‘Interferometry-by-deconvolution’ is a generalization of a 1D deconvolution method introduced by Riley and Claerbout (1976). Here we briefly review this 1D method. The 3D extension is introduced in the next section. Consider a plane wave experiment in a horizontally layered medium. At a particular depth level the total wave field is decomposed into down going and up going waves. Assuming the actual source is situated above this depth level, the total down going wave field can be seen as the illuminating wave field and the total up going wave field as its response. Subsequently, the up going wave field is deconvolved by the down going wave field. The deconvolution result is the reflection impulse response of the 1D medium below the chosen depth level. In the frequency domain, where deconvolution is replaced by division, this can be formulated as R̂+ 0 (x3,1,ω) = p̂(x3,1,ω)/ p̂(x3,1,ω), (1) where x3,1 is the x3-coordinate of the depth level at which the decomposition and division take place (in this paper the x3-axis points downwards), and p̂+ and p̂− are the down going and up going wave fields, respectively (the Fourier transform of a time-dependent function f (t) is defined as f̂ (ω) = ∫∞ −∞ f (t) exp(− jωt)dt, where j is the imaginary unit and ω denotes the angular frequency, which is taken as non-negative throughout this paper). The reflection response R̂0 (x3,1,ω) is the response that would be measured with source and receiver at x3,1 and a homogeneous half-space above x3,1. This is independent of the actual configuration above x3,1. For example, if x3,1 is chosen just below the sea-bottom, R̂0 (x3,1,ω) is the response of the medium below the sea-bottom, free of multiples related to the sea-bottom as well as to the water surface. Hence, R̂0 (x3,1,ω) obeys different boundary conditions than p̂+ and p̂−. Throughout this paper we will loosely use the term ‘deconvolution’ for division in the frequency domain (as in equation (1)). When the division is carried out for a sufficient range of frequencies, the result can be inverse Fourier transformed, yielding the time domain deconvolution result (e.g. R0 (x3,1, t)). The analogy of equation (1) with interferometry is as follows (see also Snieder, Sheiman and Calvert 2006): the righthand side is a ‘deconvolution’ of two received wave fields (instead of a correlation of two wave fields), whereas the lefthand side is the response of a virtual source at the position of a receiver (just as in interferometry). Moreover, independent of the actual source signature (transient or noise), the time domain deconvolution result R0 (x3,1, t) is an impulse response. Of course in practice the division in equation (1) should be carried out in a stabilized sense, meaning that the result becomes a band-limited impulse response. An important difference with most versions of interferometry is that equation (1) remains valid even when the medium is dissipative. Another difference is that the application of equation (1) changes the boundary conditions, as explained above. Bakulin and Calvert (2006) proposed a similar 1D deconvolution to improve their virtual source method. Snieder, Sheiman and Calvert (2006) employed a variant of this method (with source and receiver at different depth levels) to derive the impulse response of a building from earthquake data, and Mehta, Snieder and Graizer (2007) used a similar approach to estimate the near-surface properties of a dissipative medium. I N T E R F E R O M E T RYB YD E C O N V O L U T I O N : 3 D S C A L A R V E R S I O N The 1D deconvolution approach formulated by equation (1) has been extended by various authors to a multi-dimensional deconvolution method as a means for surface related and sea-bottom related multiple elimination (Wapenaar and Verschuur 1996; Ziolkowski, Taylor and Johnston 1998; Amundsen 1999; Wapenaar et al. 2000; Holvik and Amundsen 2005). In the following we derive this multi-dimensional deconvolution method along the same lines as our derivation for seismic interferometry by crosscorrelation (Wapenaar, Thorbecke and Draganov 2004). First we consider the situation for scalar fields; in the next section we generalize the derivation for vector fields. Note that when we speak of ‘fields’ we mean wave and/or diffusion fields. C © 2008 European Association of Geoscientists & Engineers, Geophysical Prospecting, 56, 419–434 Seismic and electromagnetic controlled-source interferometry 421 The starting point for our derivation is a reciprocity theorem of the convolution type for one-way scalar fields, which reads in the space-frequency domain ∫ ∂D1 { p̂A p̂B − p̂A p̂B}dx = ∫ ∂Dm { p̂A p̂B − p̂A p̂B}dx, (2) where x = (x1, x2, x3) is the Cartesian coordinate vector, ∂D1 and ∂Dm are two horizontal boundaries of infinite extent (with ∂Dm below ∂D1), and p̂+ = p̂+(x,ω) and p̂− = p̂−(x,ω) are flux-normalized down going and up going fields, respectively (see Appendix A for the derivation). The terms ‘down going’ and ‘up going’ should be interpreted in a broad sense: for diffusion fields these terms mean ‘decaying in the positive or negative x3-direction, respectively’. The subscripts A and B refer to two independent states. Equation (2) holds for lossless as well as dissipative 3D inhomogeneous media. The underlying assumptions for equation (2) are that there are no sources between ∂D1 and ∂Dm and that in the region enclosed by these boundaries the medium parameters in states A and B are identical. Above ∂D1 and below ∂Dm the medium parameters and boundary conditions in states A and B need not be the same. The condition that ∂D1 and ∂Dm are horizontal boundaries can be relaxed. Frijlink 2007 shows that under certain conditions equation (2) also holds when ∂D1 and ∂Dm are smoothly curved boundaries. Note that other variants of equation (2) exist, containing vertical derivatives of the down going and up going fields in one of the two states. This is the case, for example, when p̂+ and p̂− represent down going and up going acoustic pressure fields. Since we consider flux-normalized fields these derivatives are absent in equation (2). In the following, state B will represent the measured response of the real Earth, whereas state A will represent the new response of a redatumed source in an Earth with different boundary conditions, obtained by interferometry. Hence, state B is the actual state whereas state A is the desired state. First we discuss state B. Consider a dissipative 3D inhomogeneous Earth bounded by a free surface ∂D0, see Fig. 1(b). The source of the actual field at xS, with source spectrum ŝ(ω), is situated below ∂D0 and above the receivers. The receivers are located, for example, at the sea-bottom or in a horizontal borehole. The boundary ∂D1 is chosen an #-distance below the receivers (e.g. just below the sea-bottom) and ∂Dm is chosen below all inhomogeneities. The measured field at the receivers is represented by a 2 × 1-vector Q̂(x, xS,ω), containing for example the acoustic pressure and vertical component of the particle velocity, or the inline electric field and crossline magnetic field components. The quantities in this vector are conFigure 1. State A: the desired reflection response of the medium below ∂D1, for the situation of a non-reflecting half-space above ∂D1. State B: the actual response of the real earth, bounded by a free surface at ∂D0. The medium parameters exhibit dissipation, are 3D inhomogeneous functions of position, and below ∂D1 they are the same in both states. tinuous in the depth direction, hence, at ∂D1 (i.e. just below the receivers) we have the same Q̂(x, xS,ω). This field vector is decomposed at ∂D1 into flux-normalized down going and up going fields, according to P̂(x, xS,ω) = L̂−1Q̂(x, xS,ω), (3) where L̂−1 is a decomposition operator containing the medium parameters at ∂D1, and P̂(x, xS,ω) = ( p̂+(x, xS,ω) p̂−(x, xS,ω) )